Thermodynamics of localized magnetic moments in a Dirac conductor 
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We show that the magnetic susceptibility of a dilute ensemble of magnetic impurities in a conduc- 
tor with a relativistic electronic spectrum is non-analytic in the inverse tempertature at — >■ 0. 
We derive a general theory of this effect and construct the high-temperature expansion for the disor- 
der averaged susceptibility to any order, convergent at all tempertaures down to a possible ordering 
transition. When applied to Ising impurities on a surface of a topological insulator, the proposed 
general theory agrees with Monte Carlo simulations, and it allows us to Hnd the critical temperature 
of the ferromagnetic phase transition. 
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Collective phenomena in random ensembles of mag- 
netic impurities embedded in metallic conductors are 
caused by the long-distance exchange interaction medi- 
ated by the mobile carriers, also known as the RKKY 
exchange [Q. The RKKY interaction has a well studied 
universal structure for all metallic systems |Q, with the 
exception of the recently discovered class of two dimen- 
sional (2D) materials in which the low-energy electron 
excitations resemble massless Dirac particles: graphene 
1^, ^ , chiral metals formed at the surface of topological 
insulators |^ , and silicene Q . Recent experiments flll l 
and have reported the formation of a band gap in a 
chiral metal contaminated by magnetic impurities, point- 
ing towards magnetic ordering at the surface of the topo- 
logical insulator; also theoretical modelling suggested or- 
dering transition in some of such systems |^, |l^ . 

There are two peculiarities of the RKKY exchange in 
conductors with the Dirac-like electron spectrum which 
make it qualitatively different from usual metals: (i) the 
exchange interaction as a function of distance between 
two impurities shows the unusual decay law, (ii) the 
Friedel oscillations are either absent or comensurate with 
the lattice In the following, we shall call such in- 

teraction a Dirac-RKKY exchange. Other details of the 
RKKY exchange such as its anisotropy or whether it is 
ferromagnetic, antiferromagnetic or depolarizing may de- 
pend on the material and the symmetry of the impurity 
position in the lattice. In this Letter we propose a general 
quantitative theory of the thermodynamics of an ensem- 
ble of randomly positioned magnetic impurities interact- 
ing through the Dirac-RKKY exchange in a paramag- 
netic phase, that is at temperatures above the magnetic 
ordering temperature T^- We show that, due to a peculiar 
decay law of the Dirac-RKKY exchange, the magnetic 
susceptibility shows a strong deviation from the Curie- 
Weiss law seen as non-analiticity in its high-temperature 
expansion. 
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where cq is the Curie constant, s is the impurity spin 



Figure 1: Magnetic susceptibility of the Dirac-RKKY Ising 
model. The solid line shows the L=5 D-log Fade approximant 
calculated from the first twelve terms in the expansion (p^) 
(see Table ^ for the numerical values of the coefficients Cnj- 
The points with error bars are the results of the Monte Carlo 
simulation. The error bars show the statistical uncertainty 
arising from both MC fiuctuations and quenched disorder. 
Inset: The approximate critical temperature derived from the 
position of the singularity of the D-log Fade approximant. 
The error bars show the uncertainty due to the finite precision 
of the expansion coefficients in Table 0. The thin dotted line 
shows the best aproximation for the critical temperature. 



quantum number, Tq cx fi^/^ is a temperature scale 
dependent on the density of impurities, p, and C„ 
are numerical coefficients expressed in terms of finite- 
dimensional integrals of elementary functions, see Eqs. 
(^; (0). The expansion (^ also encodes detailed infor- 
mation about the critical point of a magnetic transition: 
the value of Tc and the susceptibility critical exponent 7 
can be extracted from the values of several coefficients 
Cn with the help of Fade approximation One ex- 

ample of a successful application of the proposed theory 
is illustrated in Fig. where the susceptibility obtained 
with the help of Eq. (|^) is compared with the Monte 
Carlo data for the archetypal Dirac-RKKY system of 
randomly positioned Ising spins with ferromagnetic I /r^ 
exchange. Numerical values of the coefficient C2 for some 
other Dirac-RKKY models are presented in Table |. 
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A pair of localized spins separated by a distance ex- 
ceeding a few lattice constants experiences the RKKY 
interaction mediated by electrons near the Fermi sur- 
face, which in Dirac conductors consists of a discrete set 
of points in the reciprocal space. In graphene there are 
two Fermi points related by time reversal symmetry [Q. 
In topological insulators where the ultrarelativistic elec- 
trons reside at the surface, there may be one Fermi point 
per face ||^. The dispersion relation of the electronic ex- 
citations near each Fermi point is e{p) = ±hvp, where p 
is the Bloch wave number of an electron relative to the 
Fermi point. Due to the discrete geometry of the Fermi 
sufrace and the linear dispersion of the excitations, the 
Dirac-RKKY interaction does not exhibit Friedel oscil- 
lations found in usual metals and decays as l/r^ as a 
function of distance r between spins ^. This decay 
law is valid in a broad range of lengths te < r < rx, 
where rE — hv/E is the length scale associated with the 
binding energy E 1 eV) of adatom, and tt = Hv/kT is 
the thermal wevelength of an electron. The most general 
form of the Dirac-RKKY Hamiltonian of a system of N 
impurities in a magnetic field is 



the spin-spin correlation function, as 



H12...N = J^-'^ 



3 hSz , 

Gij — G{Si , Sj , Ilij , , ) . 
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Here the sum is taken over all pairs of impurities ran- 
domly distributed with the average density p in the con- 



ductor plane, and 



r," I is the distance between 



a pair of adatoms with quenched positions and r^. In 
Eq. (||) we assume without loss of generality that the 
magnetic field h is coupled to the z-projection of the to- 
tal spin Sz = sf + . . . s^, and the spin of each impurity 
is assumed to be in the 2s -I- 1 dimensional representa- 
tion of SU(2). The parameter J is specific for the given 
host material and the type of impurity. Together with 
the impurity density p it defines the energy scale 



(3) 



Further material-dependent details of the RKKY interac- 
tion are encoded in the dimensionless pairwise interaction 
function Gij , which depends on two impurity spins , s ^ , 
a unit vector n^j = / Vij and two discrete quenched ran- 
dom variables ^i, ^j defining the position of each adatom 
inside the lattice unit cell [|[ |l^ (see Table | for exam- 
ples). 

The non-analiticity of the high-temperature expansion, 
Eq. (|^), results from the interplay between the peculiar 
1 /r^ Dirac-RKKY interaction and the randomness of the 
distribution of impurities in the system. The susceptibil- 
ity per spin of the system can be expressed in terms of 
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where the average is taken over both the thermal con- 
figurations of impurity spins with the Boltzmann weight 
defined by the Hamiltonian (|^), and the quenched ran- 
dom variables. The 1/r^ dependence of the Dirac-RKKY 
interaction in (^) makes it impossible to use the stan- 
dard 1/T expansion to analyse the susceptibility (Q). 
Indeed, for any temperature Tc < T < Te the ensem- 
ble contains a finite fraction of pairs, in which the spins 
are close enough to each other to be strongly correlated. 
Consider, for example, the classical ferromagnetic Ising 
model with Gy = —SiSj, where Si = ±s. Due to the 
presence of correlated pairs the high-temperature asymp- 
totics of the susceptibility splits into two contributions. 
Those spins that belong to pairs smaller than the cor- 
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strongly bound into one block spin having the length 
2s. The fraction of such spins is p - pR^ = {To/Tf/^. 
The rest of the spins are separated by distances exceed- 
ing i?o and can be regarded as an ideal gas of spins of 
length s. Then the mixture of the ideal gas of pairs and 
the ideal gas of single spins has the Curie susceptibility 
= co(l— _p)s^/T-f cop(2s)-^/2T, which deviates from the 
ideal gas susceptibility, cqs^ /T, by a non-analytic correc- 
tion 5x(xp/Toi (r)-^/3. 

A quantitative theory for the Dirac-RKKY systems 
in paramagnetic phase requires a resummation of the 
short-distance singularities appearing in the disorder av- 
erage of the observables. This is achieved by combining 
the replica method with the virial expansion of the free 
energy in the temperature-dependent gas parameter p. 
Thermodynamic properties of the system are encoded in 
the potential 



F{h, T, N; q) = -T In [Z{h, T, N)] 



(5) 



where Z{h,T, N) is the partition function of a given re- 
alization of the system of N impurities at tempearture 
T and in the presence of the magnetic field h. The inte- 
ger q is the number of identical replicas of the disordered 
system. The overline stands for the averaging over all 
quenched variables, 

- 1 f ^ 

i = l 

where A is the area of the sample, and S ~ 1 is the 
number of distinct values of the variable ^. 

The magnetic susceptibility (Q) can be written as 

x(T) = ^limlim|^|-n/i,r,7V;Q), (6) 
A* h^a q^Q oh'^ oq 
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Dirac-RKKY ensemble 


Gij 


Observable 


C2 


e/To 


Spin 1/2 Ising impurities. 






9.02 


27.1 


Spin 1/2 impurities isotropically coupled 
to the electrons in a chiral metal. 


-sfsj - (n • Si)(n • Sj) + i(n x Si)^(n x Sj)^ 




9.59 


29.7 


Sx 


1.18 


1.28 


Spin 1/2 impurities with X-Y coupling to 
the electrons in a chiral metal. 


-(n • Si)(n • Sj) + Kn X 8^)^(11 x 8^)^ 




3.17 


5.63 


Sx 


1.17 


2.16 


Spin 1/2 impurities in graphene located 
at centres of hexagons. 


-(s.-s,)cosf (C.-G), e £{0,1,2} 




-2.20 


-3.27 




3.90 


7.69 



Table I: Examples of Dirac-RKKY Hamiltonians and the corresponding values of the coefficients Ci, C2 and the generalized 
Curie temperature 0. The two dimensional host material is assumed to be in the x-y plane. The susceptibility is calculated for 
the observable given in the third column. The observable "I> in the last row defines a staggered order in which impurity spins 
in graphene are correlated with the \/3 x \/3 superlattice induced by the commensurate Friedel oscillations H FB| . 



which requires analytic continuation of the potential F 
to the positive real axis of q. In order to obtain the virial 
expansion for the susceptibility it is convenient to 
consider the grand canonical ensemble and introduce the 
thermodynamic potential 
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which is related to the potential (||) by the Legendre 
transformation, 



T = f7 + [lN, 
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The chemical potential is a g-dependent auxiliary vari- 
able, which does not have any straightforward physical 
meaning. The coefficients in SI are called the virial 
coefficients. The first three of those are: 
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where Z is a function of n quenched variables 
{r^j^ , . . . , } and {^i^ , . . . , } describing localized spins 
with indices ii, . . . ,i„ S {1, . . . ,N}, and Hi-^^, ^^i^ is the 
Hamiltonian (H) constrained to this n-subset. The trace 
is taken over all spin variables of the n-subset, and the 
symbol S represents a complete symmetrization of the 
expression over the particle indices. Note that S does 
not appear in Vi and V2 because these expressions are 
already symmetric. 

After substituting the virial expansion of SI into (H) 



one arrives at Eq. (Q), with 
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z=0 



where the integration is over n — 1 dimensionless variables 
(x2, . . . x„) = (r2, . . . , r„)/i?o- The generating function 
Q{z) in Eq. (Bh is defined by the infinite series 



Q{z) = e- 



Mi...„ = 
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zMi + ^Mi2 + ||-Mi23 + 
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Note tliat there is no integration over xi in Eq. (^), and 
tlie answer does not depend on the choice of xi by trans- 
lational invariance of the problem. Equations (|^), (||), 
and ( |lo| ) constitute the main result of this work, which 
can be applied to any random Dirac-RKKY model. 

High-temperature expansion and the critical point of 
the Ising model with random Dirac-RKKY exchange. As 
an illustration, we consider the paramagnetic suscep- 
tibihty of the ferromagnetic Ising model described by 
the pairwise interaction function Gij — —sfs^, where 
sf = ±1/2 is the Ising spin of the i-th atom. This model 
describes, for example, an ensemble of easy-axis magnetic 
impurities on the sufrace of a topological insulator [^Sj; 
it also describes structural phase transitions in graphene 
1^. The only random quenched variables in this model 
are the positions of impurities in the sample. The first 
two virial coefficients found from Eqs. ffl), dlG) are 



Ci = 1 



C2 



1 

2./0 



2T:xdx = 9.01845 . 



. 1 + e ^ 

Similarly, all higher-order C„ are expressed as integrals 
of elementary functions of increasing complexity. The 
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Ci 1 Cs 3.128(2) X 10" C9 5.612(4) x 10*^ 

C2 9.0184534. . . Cb 2.103(1) x 10"* do 3.328(3) x lO'^ 

C3 6.578(2) X 10^ C7 1.3937(5) x 10^ di 1.857(3) x 10*^ 

C4 4.582(1) X 10^ Cs 9.007(4) x 10^ C12 1.048(7) x 10'^ 



Table II: Numerical values of the coefficients in the high- 
temperature expansion (jl|) for the susceptibility of the fer- 
romagnetic Ising model. 



numerical values of the first twelve coefficients C„ are 
given in Table ||. 

At low enougli temperature the Ising model with a fer- 
romagnetic exchange undergoes a phase transition into a 
ferromagnetically ordered state. In order to extract the 
properties of the critical point from the high-temperature 
data in Table |l| we use the D-log Pade technique 
to investigate the function f{z) = J2n>i CnZ^~^ , where 
Cn are the same coefficients as appear in Eq. (^, and 
z = (To/T)^/'^. More precisely, we construct the first 
few [L/L + 1] Pade approximants, < L < 5, to the 
function g{z) — f'{z)/f(z) and analyse the structure 
of the singularity of the approximants in the vicinity 
of the hypothetical critical point. For the L — h ap- 
proximant, which requires the knowledge of the first 
twelve coefficients C„, we solve the differential equation 
d\x\\f {z)] / dz = [L/L + l]g(^) with the initial condition 
/(O) = and compare the result with the susceptibil- 
ity obtained from Monte Carlo simulations on an ensem- 
ble of 10^ impurities averaged over 100 configurations of 
quenched disorder, see Fig. 0. The simulations are per- 
formed with the classical worm algorithm |0 using the 
ALPS libraries ||l8|. Since the numerical values of C„ 
have been calculated with finite precision, the position of 
the singularity of each approximant is found with some 
uncertainty (see the inset of Fig. (Q)). For the L — 5 ap- 
proximant we find Tc/Tq — 14.5(5), which is close to the 
value found from numerical simulations in ^ . We 

also extract the value of the susceptibility critical expo- 
nent 7 by calculating the residue of the L — 5 Pade ap- 
proximant at the critical point and find that 7 — 1.4(2). 

We conclude this report with the discussion of several 
other examples of the Dirac-RKKY systems identified in 
earlier literature. A detailed analysis of criticality in such 
systems is beyond the scope of this letter, however there 
is enough interesting information already contained in 
the leading coefficients Ci and C2 of the expansion (|^) . 
Not only do these coefficient describe an experimentally 
measurable deviation of xC^) from the Curie law, but 
also they can be used to extract a generalized Curie tem- 
perature Q — sign(Ci) x (|C2|/Ci)^/^To, which gives an 
estimate of Tc if the susceptibility is calculated for the 
observable corresponding to the order parameter. For 
spin systems having competing orders one can compare 
the generalized Curie temperatures extracted from the 



susceptibilities for the suspected order parameters for 
which ordering is more likely to occur. Such invorma- 
tion is given in Table |. The first two rows describe two 
limiting cases of a more general model derived in [T^ 
for the RKKY exchange in chiral metals. The authors 
of ||l^ observe that a strongly anisotropic in-plane ex- 
change tends to destroy the magnetic order enforced by 
the out-of-plane ferromagnetic coupling. They also con- 
jecture that as a function of the anisotropy parameter 
the system undergoes a quantum phase transition into a 
spin-glass state. We do not find any signature of such a 
transition in the generalized Curie temperature for both 
in- and out-of-plane susceptibility. 

The third row of Table | describes a Dirac-RKKY 
model of magnetic adatoms at the centers of carbon 
hexagons in graphene |p^ . The quenched parameter f 
labels the three distinct positions of the atom in the 
\/3x -\/3 superlattice formed by the commensurate Friedel 
oscillations of the RKKY interaction. In this case ferro- 
magnetic ordering is replaced by a staggered state with 
the order parameter 

$ = ^[s^ cos(27r^,/3) -I- sf sin(27r^,/3)]. 

i 

Together with the detailed analysis of the Ising system 
these examples demonstrate how the proposed generic 
high-temperature expansion, Eqs.(||), (^ and ( |lo|) can 
be used to describe the magnetic properties of disordered 
Dirac-RKKY systems. 
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